#ifndef SLEPBASE_H
#define SLEPBASE_H

#ifdef __cplusplus
extern "C" {
#endif

#include <math.h>

#define SLEPMAX(x,y) (((x)>(y))?(x):(y))
#define SLEPMIN(x,y) (((x)<(y))?(x):(y))
#define SLEPABS(x) ((x)<0.0?-(x):(x))
#define SLEPDoubleApproxEquals(x,y) (SLEPABS((x)-(y))<0.0001*SLEPMAX((x),(y))?1:0)

typedef enum{
	SlepLeft,
	SlepRight
}slep_side;

typedef enum{
	SlepCompressedSparseColumn,
	SlepCompressedSparseRow
}slep_sparse_format;

typedef enum{
	SlepSparse,
	SlepDense,
	SlepSymmetric,
	SlepTriangular,
	SlepBanded
}slep_sparsity;

typedef enum{
	SlepTranspose,
	SlepNoTranspose
}slep_transposity;

typedef enum{
	SlepUpper,
	SlepLower,
	SlepBandedUpper,
	SlepBandedLower,
	SlepBandedBoth,
	SlepFull
}slep_uplobanded;

typedef struct{
	double *values;
	int *ColumnsORRows;
	int *Index;
	slep_sparse_format format;
}sparse_matrix_base; // Zero Based Indexing Compressed-Sparse-Row Format

typedef struct{
	union{
		sparse_matrix_base spmtrx;
		double *mtrx;
	}data;

	/* This has value in {SlepDense, SlepSparse, SlepSymmetric, SlepTriangular} */
	slep_sparsity sparsity;

	/* This has value in {SlepUpper, SlepLower, SlepBandedUpper, SlepBandedLower, SlepBandedBoth, SlepFull} */
	slep_uplobanded uplobanded;

	/* M Rows */
	int m;

	/* N Columns */
	int n;

	/* Leading Dimension (used for Dense Matrices and Symmetric or Banded Matrices) */
	int ld;

	/* Bandwidth */
	int kl;
	int ku;

}slep_matrix;

typedef struct{
	double *values;
	int *index;
	int nnz;
}sparse_vector_base;

typedef struct{
	union{
		sparse_vector_base spvctr;
		double *vctr;
	}data;

	/* This has value in {SlepDense, SlepSparse} */
	slep_sparsity sparsity;

	int inc;
	int dim;
}slep_vector;

slep_matrix* createCSCSlepMatrix(int m, int n, double *values, int *rows, int *columnIndex);
slep_matrix* createCSRSlepMatrix(int m, int n, double *values, int *columns, int *rowIndex);
slep_matrix* createDenseSlepMatrix(int m, int ld, int n, double *values);

slep_matrix* createSymmetricSlepMatrix(int n, double *values, slep_uplobanded uplobanded, int bandwidth);
slep_matrix* createTriangularSlepMatrix(int n, double *values, slep_uplobanded uplobanded, int bandwidth);

slep_vector* createSparseSlepVector(int nnz, int *index, double *values);
slep_vector* createDenseSlepVector(int dimension, double *values, int inc);

int slepFreeVector(slep_vector* vctr);
slep_vector* slepCopyToNewVector(slep_vector* vctr);

int slepFreeMatrix(slep_matrix* mtrx);
int findNnz(slep_matrix* mtrx);

/* Should Implement these methods */
slep_vector* slepAllocDenseVector(int size);
slep_vector* slepAllocSparseVector(int nnz);
slep_matrix* slepAllocDenseMatrix(int m,int n);
slep_matrix* slepAllocSparseMatrix(int m,int nnz);
slep_vector* slepAllocVector(int memory, slep_sparsity sparsity);


/*--------------------------- Copied from Old Slep.h File ---------------------------*/

typedef struct{
	double* x;
	double* funVal;
	double* ValueL;
	char type[20]; // strcpy the type here
	double c; // Used for logistic regression only
	int totIter; // Total Iterations that happened in this call
	int errorCode; // 0 for no error
}FUNVAL;

typedef enum{
	SlepOPT_x0 = 0x1,
	SlepOPT_c0 = 0x2,
	SlepOPT_maxIter = 0x4,
	SlepOPT_tol = 0x8,
	SlepOPT_mu = 0x10,
	SlepOPT_nu = 0x20,
	SlepOPT_rsL2 = 0x40,
	SlepOPT_indices = 0x80,
	SlepOPT_sWeight = 0x100,
	SlepOPT_X0 = 0x200,
	SlepOPT_q = 0x400,
	SlepOPT_gWeights = 0x800,
	SlepOPT_fusedPenalty = 0x1000,
	SlepOPT_rStartNum = 0x2000
}OPTS_REGISTRY;

typedef struct{

	OPTS_REGISTRY REGISTER;

	//Starting points
	slep_vector* x0; // Starting x
	slep_matrix* X0; // Starting X (matrix)
	slep_vector* c0; // Used in Logistic regression only
	double *mt_c0; // Used in mtLogistic regression only - c0 not used here
	double *mc_c0; // Used in mcLogistic regression only - c0 not used here

	int init; // Init takes values 0,1,2 default:2

	// Termination
	int maxIter;
	int rStartNum;
	double tol;
	int tFlag;

	//Normalization
	int nFlag;
	slep_vector* mu;
	slep_vector* nu;

	//Method
	int method;
	int mFlag;
	int lFlag;

	//Regularization
	int rFlag;
	double rsL2;
	double fusedPenalty;

	//Group and others
	int* ind;
	int ind_size;
	double q;
	slep_vector *gWeights;

	double sWeight[2]; // registry set to 'sWeight'
	int mc_k;
	// 2 more to come

}OPTS;

#define setREGISTER(opts,field) ((opts)->REGISTER=(OPTS_REGISTRY)(((opts)->REGISTER)|(field)))
#define unsetREGISTER(opts,field) ((opts)->REGISTER=(OPTS_REGISTRY)(((opts)->REGISTER)&(~(field))))

void initOPTS(OPTS *opt);

//int hasGot(OPTS *opt,OPTS_REGISTRY field);
#define hasGot(opts,field) ((((opts)->REGISTER)&(field))==(field))

//void setOPTx0(OPTS *opt,slep_vector* x0);
#define setOPTx0(opts,x0_v)	{\
								if(!hasGot(opts,SlepOPT_x0))\
									((opts)->REGISTER) |= SlepOPT_x0;\
								((opts)->x0) = (x0_v);\
							}

//void setOPTX0(OPTS *opt,slep_matrix* X0);
#define setOPTX0(opts,X0_v)	{\
								if(!hasGot(opts,SlepOPT_X0))\
									((opts)->REGISTER) |= SlepOPT_X0;\
								((opts)->X0) = (X0_v);\
							}


//void setOPTc0(OPTS *opt,double c0);
#define setOPTc0(opts,c0_v)	{\
								if(!hasGot(opts,SlepOPT_c0))\
									((opts)->REGISTER) |= SlepOPT_c0;\
								((opts)->c0) = (c0_v);\
							}

#define setOPTmu(opts,mu_v)	{\
								if(!hasGot(opts,SlepOPT_mu))\
									((opts)->REGISTER) |= SlepOPT_mu;\
								((opts)->mu) = (mu_v);\
							}


#define setOPTnu(opts,nu_v)	{\
								if(!hasGot(opts,SlepOPT_nu))\
									((opts)->REGISTER) |= SlepOPT_nu;\
								((opts)->nu) = (nu_v);\
							}

#define setOPTsWeight(opts,sWeight_v) {\
										if(!hasGot(opts,SlepOPT_sWeight))\
											((opts)->REGISTER) |= SlepOPT_sWeight;\
										((opts)->sWeight) = (sWeight_v);\
									}

#define setOPTind(opts,ind_v,ind_size_v)	{\
										if(!hasGot(opts,SlepOPT_indices))\
												((opts)->REGISTER) |= SlepOPT_indices;\
											((opts)->ind) = (ind_v);\
											((opts)->ind_size) = (ind_size_v);\
										}

#define setOPTq(opts,q_v)		{\
										if(!hasGot(opts,SlepOPT_q))\
											((opts)->REGISTER) |= SlepOPT_q;\
										((opts)->q) = (q_v);\
									}

#define setOPTgWeights(opts,sWeight_v) {\
										if(!hasGot(opts,SlepOPT_gWeights))\
											((opts)->REGISTER) |= SlepOPT_gWeights;\
										((opts)->gWeights) = (gWeights_v);\
									}


void sll_opts(OPTS* opts);
void initNormalization(slep_matrix* A, OPTS* opts);


/*----------END----------------- Copied from Old Slep.h File ------------END---------------*/


/* =============================== SLEP BLAS ========================================= */
/* Level 1 BLAS */

/* [dest = orig] */
void slep_dvcopy(slep_vector* orig, slep_vector* dest);
void slep_dmcopy(slep_matrix* orig, slep_matrix* dest);
slep_matrix* slepCopyToNewMatrix(slep_matrix* mtrx);
slep_vector* slepCopyToNewVector(slep_vector* vctr);

/* [x = alpha * x] */
void slep_dscal(double alpha, slep_vector* x);

/* [y = y + alpha * x] */
void slep_daxpy(double alpha, slep_vector* x, slep_vector* y);

/* [y = alpha * x + beta * y] */
void slep_daxpby(double alpha, slep_vector* x, double beta, slep_vector* y);

/* [= x.y] */
double slep_ddot(slep_vector* x, slep_vector* y);

/* [=sqrt(x.x)] */
double slep_dnrm2(slep_vector* x);

/* [=|x|_p (p-norm of x)] */
double slep_dpnrm(slep_vector* x, double p);

/* [=|x| (1-Norm)] */
double slep_dasum(slep_vector* x);

/* [=sum(x_i) (Sum of all elements of x)] */
double slep_dsum(slep_vector* x);

/* [=i such that |x_i| is maximum ] */
int slep_idamax(slep_vector* x);

/* [=|x_i| such that |x_i| is maximum ] */
double slep_damax(slep_vector* x);

/* [= x_i (such that x_i is maximum) ] */
double slep_dmax(slep_vector* x);


/* [=||A|| (Frobenius-Norm)] */
double slep_dfrbnrm(slep_matrix* a);

/* [= trace ( A` . B` ) [ ` = transposity] */
double slep_dtraceAB(slep_transposity transA, slep_matrix* a, slep_transposity transB, slep_matrix* b);

/* [y = x.diag(y) (or diag(x).y as both are same)] */
void slep_xdiagy(slep_vector* x, slep_vector* y);


/* [y = {diag(x)^-1}.y ] */
void slep_invdiagxy(slep_vector* x, slep_vector* y);





/* Level 2 BLAS */

/* [y = alpha.A`.x + beta.y] */
void slep_daAxpby(slep_transposity trans, double alpha,
				  slep_matrix* A, slep_vector* x,
				  double beta, slep_vector* y);


/* [x = A`.x (` implies SlepTranspose or SlepNoTranspose)]*/
void slep_dtrmv(slep_transposity trans,	slep_matrix* A,	slep_vector* x);

/* [x = {A`^(-1)} . x (` implies SlepTranspose or SlepNoTranspose)]*/
void slep_dtrsv(slep_transposity trans,	slep_matrix* A,	slep_vector* x);


void slep_dger(double alpha, slep_vector* x, slep_vector* y,
			   slep_matrix* A);

/* Level 3 BLAS */

void slep_dgemm(slep_transposity transA, slep_transposity transB, double alpha,
				slep_matrix* A, slep_matrix* B, double beta, slep_matrix* C);


void slep_dsymm(double alpha, slep_matrix* A, slep_matrix* B,
				double beta, slep_matrix* C);

void slep_dtrmm(slep_side side, slep_transposity transA, double alpha,
				slep_matrix* A, slep_matrix* B);


void slep_dtrsm(slep_side side, slep_transposity transA, double alpha,
				slep_matrix* A, slep_matrix* B);



/* Other useful functions */

void slep_zerom(slep_matrix *mtrx);
void slep_zerov(slep_vector *vctr);


/* =============================== SLEP BLAS ========================================= */

/* =============================== SLEP UTIL ========================================= */
#ifdef WIN32
#include <windows.h>

typedef struct {
    LARGE_INTEGER start;
    LARGE_INTEGER stop;
} stopWatch;

void startTimer( stopWatch *timer);
void stopTimer( stopWatch *timer);
double LIToSecs( LARGE_INTEGER * L);
double getElapsedTime( stopWatch *timer);
#endif



double Random(void);
void   PlantSeeds(long x);
void   GetSeed(long *x);
void   PutSeed(long x);
void   SelectStream(int index);
void   TestRandom(void);

long Bernoulli(double p);
long Binomial(long n, double p);
long Equilikely(long a, long b);
long Geometric(double p);
long Pascal(long n, double p);
long Poisson(double m);

double Uniform(double a, double b);
double Exponential(double m);
double Erlang(long n, double b);
double Normal(double m, double s);
double Lognormal(double a, double b);
double Chisquare(long n);
double Student(long n);


void slep_randn(double* x, int n, double mu, double sigma);


/* =============================== SLEP UTIL ========================================= */


/* ============================= SLEP FILE & IO ====================================== */

#define SLEPWD "C:\\Work\\testfolder\\"

slep_vector* slepReadDenseVector(char *filename);
slep_matrix* slepReadDenseMatrix(char *filename);
void slepWriteDenseVector(char* filename, slep_vector* vctr);
void slepVerboseWriteDenseVector(char* filename, slep_vector* vctr);

void slepWriteDenseMatrix(char* filename, slep_matrix* mtrx);
void slepVerboseWriteDenseMatrix(char* filename, slep_matrix* mtrx);

void slepVerbConsoleWrDMat(slep_matrix* mtrx);
void slepVerbConsoleWrDVec(slep_vector* vctr);


int LeastRTestImplementation(slep_matrix* A, slep_vector* y, double z, OPTS *opts,
			slep_vector* x, slep_vector* funVal, slep_vector* valueL);

/* ============================= SLEP FILE & IO ====================================== */




/* ============================= SLEP MACROS ====================================== */

#include <assert.h>


#define SLEP_DSCAL(alpha,x,i) 	{\
							if ( (alpha)==1.0)\
								;\
							else if ( (alpha)==0. &&  (x)->sparsity == SlepDense &&  (x)->dim == 1)\
								for((i)=0; (i) < ( (x)->dim) ; ++(i))\
									( (x)->data.vctr)[(i)] = 0.0;\
							else if ( (alpha)==0. &&  (x)->sparsity == SlepSparse )\
								for((i)=0; (i) < ( (x)->data.spvctr.nnz) ; ++(i))\
									( (x)->data.spvctr.values)[(i)] = 0.0;\
							else if( (x)->sparsity == SlepSparse)\
								for((i)=0; (i) < ( (x)->data.spvctr.nnz) ; (i)++)\
									( (x)->data.spvctr.values[(i)]) *=  (alpha);\
							else if(x->sparsity == SlepDense)\
								for((i)=0; (i) < ( (x)->dim) ; (i)++)\
									( (x)->data.vctr[(i)*( (x)->inc)]) *=  (alpha);\
						}



#define DMEMCPY(dest, orig,size,i)	{\
										for( (i)=0; (i) < (size); ++(i) )\
											(*((dest)+(i))) = (*((orig)+(i)));\
									}



#define DMEMSET(x,value,size,i)		{\
										for((i)=0; (i) < (size) ; ++(i))\
											(*((x)+(i))) = (value);\
									}


#define MVDIMCHECK(A,x)             {\
										if((A)->sparsity== SlepDense && (x)->sparsity==SlepDense)\
											assert((A)->n == (x)->dim);\
										else if( (A)->sparsity == SlepSparse && (x)->sparsity== SlepDense)\
											assert( (A)->n == (x)->dim);\
										else if( (A)->sparsity == SlepSparse && (x)->sparsity== SlepSparse)\
											assert( (A)->n >= ((x)->data.spvctr.index[(x)->data.spvctr.nnz -1]+1) );\
										else if ((A)->sparsity== SlepDense && (x)->sparsity== SlepSparse)\
											assert( (A)->n >= ((x)->data.spvctr.index[(x)->data.spvctr.nnz -1]+1) );\
   									 }


#define NORMALIZEDMTV(A,b,ATb,opts,b1)		{\
											if( (opts)->nFlag==0)\
												slep_daAxpby(SlepTranspose,1.0,(A),(b),0.,(ATb));\
											else if( (opts)->nFlag == 1){\
												slep_daAxpby(SlepTranspose,1.0,(A),(b),-(slep_dsum((b))),(opts)->mu);\
												slep_invdiagxy( (opts)->nu,(ATb));\
											}else{\
												slep_dvcopy((b),(b1) );\
												slep_invdiagxy( (opts)->nu, (b1));\
												slep_dvcopy((opts)->mu,(ATb) );\
												slep_daAxpby(SlepTranspose,1.,(A),(b1),-(slep_dsum((b1))),(ATb));\
											}	\
									}



#define NORMALIZEDMV(A,x,Ax,opts,x1,i)	{\
											if( opts->nFlag == 0)\
												slep_daAxpby(SlepNoTranspose, 1., (A) , (x) , 0., (Ax) );\
											else if( (opts)->nFlag == 1){\
												slep_dvcopy( (x) , (x1) );\
												slep_invdiagxy((opts)->nu,(x1));\
												assert( (Ax)->sparsity == SlepDense);\
												DMEMSET( (Ax)->data.vctr,1.0,(Ax)->dim, (i) );\
												slep_daAxpby(SlepNoTranspose, 1.0, (A) , (x1),  slep_ddot( (opts)->mu, (x1) ), (Ax) );\
											}else{\
												DMEMSET( (Ax)->data.vctr,1.0,(Ax)->dim, (i) );\
												slep_daAxpby(SlepNoTranspose, 1.0, (A) , (x), slep_ddot((opts)->mu, (x) ), (Ax) );\
												slep_invdiagxy((opts)->nu, (Ax) ); \
											}\
										}



/* ============================= SLEP MACROS ====================================== */


/*----------------------------- slepblaslevel1.1 ---------------------------- */
void slep_dmscal(double alpha, slep_matrix* X);
void slep_daXpbY (slep_transposity transX, double alpha, slep_matrix* X, double beta, slep_matrix* Y);
typedef enum{
	SlepNoError=0,
	SlepDimensionMismatchError=1,
	SlepNotSupportedYet=2
}slep_error;
/*----------------------------- slepblaslevel1.1 ---------------------------- */

#ifdef __cplusplus
}
#endif

#endif //SLEPBASE_H